Harvest Index a key trait to select tolerant potato genotypes (Solanum tuberosum) under drought stress condition

Journal of Crop Science and Biotechnology

1 Setup

source("files/setup.r")
library(gtsummary)
library(factoextra)
library(corrplot)

devtools::session_info()

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.2.2 (2022-10-31 ucrt)
 os       Windows 10 x64 (build 22621)
 system   x86_64, mingw32
 ui       RTerm
 language (EN)
 collate  English_United States.utf8
 ctype    English_United States.utf8
 tz       America/Buenos_Aires
 date     2023-02-04
 pandoc   2.19.2 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package       * version  date (UTC) lib source
 agricolae       1.3-5    2021-06-06 [1] CRAN (R 4.2.2)
 AlgDesign       1.2.1    2022-05-25 [1] CRAN (R 4.2.0)
 askpass         1.1      2019-01-13 [1] CRAN (R 4.2.2)
 assertthat      0.2.1    2019-03-21 [1] CRAN (R 4.2.2)
 backports       1.4.1    2021-12-13 [1] CRAN (R 4.2.0)
 boot            1.3-28   2021-05-03 [2] CRAN (R 4.2.2)
 broom           1.0.2    2022-12-15 [1] CRAN (R 4.2.2)
 broom.helpers   1.11.0   2023-01-06 [1] CRAN (R 4.2.2)
 cachem          1.0.6    2021-08-19 [1] CRAN (R 4.2.2)
 callr           3.7.3    2022-11-02 [1] CRAN (R 4.2.2)
 cellranger      1.1.0    2016-07-27 [1] CRAN (R 4.2.2)
 cli             3.6.0    2023-01-09 [1] CRAN (R 4.2.2)
 cluster         2.1.4    2022-08-22 [2] CRAN (R 4.2.2)
 codetools       0.2-18   2020-11-04 [2] CRAN (R 4.2.2)
 colorspace      2.1-0    2023-01-23 [1] CRAN (R 4.2.2)
 combinat        0.0-8    2012-10-29 [1] CRAN (R 4.2.0)
 corrplot      * 0.92     2021-11-18 [1] CRAN (R 4.2.2)
 cowplot       * 1.1.1    2020-12-30 [1] CRAN (R 4.2.2)
 crayon          1.5.2    2022-09-29 [1] CRAN (R 4.2.2)
 curl            5.0.0    2023-01-12 [1] CRAN (R 4.2.2)
 DBI             1.1.3    2022-06-18 [1] CRAN (R 4.2.2)
 dbplyr          2.3.0    2023-01-16 [1] CRAN (R 4.2.2)
 devtools      * 2.4.5    2022-10-11 [1] CRAN (R 4.2.2)
 digest          0.6.31   2022-12-11 [1] CRAN (R 4.2.2)
 dplyr         * 1.0.10   2022-09-01 [1] CRAN (R 4.2.2)
 DT              0.27     2023-01-17 [1] CRAN (R 4.2.2)
 ellipsis        0.3.2    2021-04-29 [1] CRAN (R 4.2.2)
 emmeans         1.8.4-1  2023-01-17 [1] CRAN (R 4.2.2)
 estimability    1.4.1    2022-08-05 [1] CRAN (R 4.2.1)
 evaluate        0.20     2023-01-17 [1] CRAN (R 4.2.2)
 factoextra    * 1.0.7    2020-04-01 [1] CRAN (R 4.2.2)
 FactoMineR    * 2.7      2022-12-14 [1] CRAN (R 4.2.2)
 fansi           1.0.4    2023-01-22 [1] CRAN (R 4.2.2)
 fastmap         1.1.0    2021-01-25 [1] CRAN (R 4.2.2)
 flashClust      1.01-2   2012-08-21 [1] CRAN (R 4.2.0)
 forcats       * 0.5.2    2022-08-19 [1] CRAN (R 4.2.2)
 fs              1.6.0    2023-01-23 [1] CRAN (R 4.2.2)
 gargle          1.2.1    2022-09-08 [1] CRAN (R 4.2.2)
 generics        0.1.3    2022-07-05 [1] CRAN (R 4.2.2)
 ggplot2       * 3.4.0    2022-11-04 [1] CRAN (R 4.2.2)
 ggrepel         0.9.2    2022-11-06 [1] CRAN (R 4.2.2)
 glue            1.6.2    2022-02-24 [1] CRAN (R 4.2.2)
 googledrive   * 2.0.0    2021-07-08 [1] CRAN (R 4.2.2)
 googlesheets4 * 1.0.1    2022-08-13 [1] CRAN (R 4.2.2)
 gsheet        * 0.4.5    2020-04-07 [1] CRAN (R 4.2.2)
 gt              0.8.0    2022-11-16 [1] CRAN (R 4.2.2)
 gtable          0.3.1    2022-09-01 [1] CRAN (R 4.2.2)
 gtsummary     * 1.7.0    2023-01-13 [1] CRAN (R 4.2.2)
 haven           2.5.1    2022-08-22 [1] CRAN (R 4.2.2)
 highr           0.10     2022-12-22 [1] CRAN (R 4.2.2)
 hms             1.1.2    2022-08-19 [1] CRAN (R 4.2.2)
 htmltools       0.5.4    2022-12-07 [1] CRAN (R 4.2.2)
 htmlwidgets     1.6.1    2023-01-07 [1] CRAN (R 4.2.2)
 httpuv          1.6.8    2023-01-12 [1] CRAN (R 4.2.2)
 httr            1.4.4    2022-08-17 [1] CRAN (R 4.2.2)
 huito         * 0.2.2    2023-01-24 [1] local
 inti          * 0.6.0    2023-01-24 [1] CRAN (R 4.2.2)
 jsonlite        1.8.4    2022-12-06 [1] CRAN (R 4.2.2)
 klaR            1.7-1    2022-06-27 [1] CRAN (R 4.2.2)
 knitr         * 1.41     2022-11-18 [1] CRAN (R 4.2.2)
 labelled        2.10.0   2022-09-14 [1] CRAN (R 4.2.2)
 later           1.3.0    2021-08-18 [1] CRAN (R 4.2.2)
 lattice         0.20-45  2021-09-22 [2] CRAN (R 4.2.2)
 leaps           3.1      2020-01-16 [1] CRAN (R 4.2.2)
 lifecycle       1.0.3    2022-10-07 [1] CRAN (R 4.2.2)
 lme4            1.1-31   2022-11-01 [1] CRAN (R 4.2.2)
 lubridate       1.9.1    2023-01-24 [1] CRAN (R 4.2.2)
 magick        * 2.7.3    2021-08-18 [1] CRAN (R 4.2.2)
 magrittr        2.0.3    2022-03-30 [1] CRAN (R 4.2.2)
 MASS            7.3-58.1 2022-08-03 [2] CRAN (R 4.2.2)
 Matrix          1.5-1    2022-09-13 [2] CRAN (R 4.2.2)
 memoise         2.0.1    2021-11-26 [1] CRAN (R 4.2.2)
 mime            0.12     2021-09-28 [1] CRAN (R 4.2.0)
 miniUI          0.1.1.1  2018-05-18 [1] CRAN (R 4.2.2)
 minqa           1.2.5    2022-10-19 [1] CRAN (R 4.2.2)
 mnormt          2.1.1    2022-09-26 [1] CRAN (R 4.2.1)
 modelr          0.1.10   2022-11-11 [1] CRAN (R 4.2.2)
 multcomp        1.4-20   2022-08-07 [1] CRAN (R 4.2.2)
 multcompView    0.1-8    2019-12-19 [1] CRAN (R 4.2.2)
 munsell         0.5.0    2018-06-12 [1] CRAN (R 4.2.2)
 mvtnorm         1.1-3    2021-10-08 [1] CRAN (R 4.2.0)
 nlme            3.1-160  2022-10-10 [2] CRAN (R 4.2.2)
 nloptr          2.0.3    2022-05-26 [1] CRAN (R 4.2.2)
 openssl         2.0.5    2022-12-06 [1] CRAN (R 4.2.2)
 pillar          1.8.1    2022-08-19 [1] CRAN (R 4.2.2)
 pkgbuild        1.4.0    2022-11-27 [1] CRAN (R 4.2.2)
 pkgconfig       2.0.3    2019-09-22 [1] CRAN (R 4.2.2)
 pkgload         1.3.2    2022-11-16 [1] CRAN (R 4.2.2)
 prettyunits     1.1.1    2020-01-24 [1] CRAN (R 4.2.2)
 processx        3.8.0    2022-10-26 [1] CRAN (R 4.2.2)
 profvis         0.3.7    2020-11-02 [1] CRAN (R 4.2.2)
 promises        1.2.0.1  2021-02-11 [1] CRAN (R 4.2.2)
 ps              1.7.2    2022-10-26 [1] CRAN (R 4.2.2)
 psych         * 2.2.9    2022-09-29 [1] CRAN (R 4.2.2)
 purrr         * 1.0.1    2023-01-10 [1] CRAN (R 4.2.2)
 questionr       0.7.7    2022-01-31 [1] CRAN (R 4.2.2)
 R6              2.5.1    2021-08-19 [1] CRAN (R 4.2.2)
 rappdirs        0.3.3    2021-01-31 [1] CRAN (R 4.2.2)
 Rcpp            1.0.10   2023-01-22 [1] CRAN (R 4.2.2)
 readr         * 2.1.3    2022-10-01 [1] CRAN (R 4.2.2)
 readxl          1.4.1    2022-08-17 [1] CRAN (R 4.2.2)
 remotes         2.4.2    2021-11-30 [1] CRAN (R 4.2.2)
 reprex          2.0.2    2022-08-17 [1] CRAN (R 4.2.2)
 rlang           1.0.6    2022-09-24 [1] CRAN (R 4.2.2)
 rmarkdown       2.20     2023-01-19 [1] CRAN (R 4.2.2)
 rstudioapi      0.14     2022-08-22 [1] CRAN (R 4.2.2)
 rvest           1.0.3    2022-08-19 [1] CRAN (R 4.2.2)
 sandwich        3.0-2    2022-06-15 [1] CRAN (R 4.2.2)
 scales          1.2.1    2022-08-20 [1] CRAN (R 4.2.2)
 scatterplot3d   0.3-42   2022-09-08 [1] CRAN (R 4.2.1)
 sessioninfo     1.2.2    2021-12-06 [1] CRAN (R 4.2.2)
 shiny         * 1.7.4    2022-12-15 [1] CRAN (R 4.2.2)
 showtext        0.9-5    2022-02-09 [1] CRAN (R 4.2.2)
 showtextdb      3.0      2020-06-04 [1] CRAN (R 4.2.2)
 stringi         1.7.12   2023-01-11 [1] CRAN (R 4.2.2)
 stringr       * 1.5.0    2022-12-02 [1] CRAN (R 4.2.2)
 survival        3.4-0    2022-08-09 [2] CRAN (R 4.2.2)
 sysfonts        0.8.8    2022-03-13 [1] CRAN (R 4.2.2)
 TH.data         1.1-1    2022-04-26 [1] CRAN (R 4.2.2)
 tibble        * 3.1.8    2022-07-22 [1] CRAN (R 4.2.2)
 tidyr         * 1.3.0    2023-01-24 [1] CRAN (R 4.2.2)
 tidyselect      1.2.0    2022-10-10 [1] CRAN (R 4.2.2)
 tidyverse     * 1.3.2    2022-07-18 [1] CRAN (R 4.2.2)
 timechange      0.2.0    2023-01-11 [1] CRAN (R 4.2.2)
 tzdb            0.3.0    2022-03-28 [1] CRAN (R 4.2.2)
 urlchecker      1.0.1    2021-11-30 [1] CRAN (R 4.2.2)
 usethis       * 2.1.6    2022-05-25 [1] CRAN (R 4.2.2)
 utf8            1.2.2    2021-07-24 [1] CRAN (R 4.2.2)
 vctrs           0.5.2    2023-01-23 [1] CRAN (R 4.2.2)
 withr           2.5.0    2022-03-03 [1] CRAN (R 4.2.2)
 xfun            0.36     2022-12-21 [1] CRAN (R 4.2.2)
 xml2            1.3.3    2021-11-30 [1] CRAN (R 4.2.2)
 xtable          1.8-4    2019-04-21 [1] CRAN (R 4.2.2)
 yaml            2.3.7    2023-01-23 [1] CRAN (R 4.2.2)
 zoo             1.8-11   2022-09-17 [1] CRAN (R 4.2.2)

 [1] C:/Users/floza/AppData/Local/R/win-library/4.2
 [2] C:/Program Files/R/R-4.2.2/library

──────────────────────────────────────────────────────────────────────────────

2 Googlesheets connect

url <- "https://docs.google.com/spreadsheets/d/1dfgpmCKdPmxRHozrZp0iE_xMGsvKTIcztDpMWYSEGaY"
# browseURL(url)
gs <- as_sheets_id(url)

3 Trait rename

geno <- c("CIP720088" = "G01"
          , "CIP392797.22" = "G02"
          , "CIP397077.16" = "G03"
          , "CIP398192.213" = "G04"
          , "CIP398180.612" = "G05"
          , "CIP398208.704" = "G06"
          , "CIP398098.119" = "G07"
          , "CIP398190.89" = "G08"
          , "CIP398192.592" = "G09"
          , "CIP398201.510" = "G10"
          , "CIP398203.244" = "G11"
          , "CIP398203.5" = "G12"
          , "CIP398208.219" = "G13"
          , "CIP398208.33" = "G14"
          , "CIP398208.620" = "G15")

traits <- c("SPAD_29" = "Chlorophyll concentration (SPAD) at 29 dap"
            , "SPAD_59" = "Chlorophyll concentration (SPAD) at 59 dap"
            , "SPAD_76" = "Chlorophyll concentration (SPAD) at 76 dap"
            , "SPAD_83" = "Chlorophyll concentration (SPAD) at 83 dap"
            , "HGT" = "Plant height (cm)"
            , "RWC" = "Relative water content (%)"
            , "LOP" = "Leaf osmotic potential (MPa)"
            , "LDW" = "Leaf dry weight (g)"
            , "SDW" = "Stem dry weight (g)"
            , "RDW" = "Root dry weight (g)"
            , "TDW" = "Tuber dry weight (g)"
            , "NTUB" = "Tuber number (N°)" 
            , "TRS" = "Total transpiration (mL)"
            , "LFA" = "Leaf area (cm^2^)"
            , "RTL" = "Root length (cm)"
            , "TDB" = "Total dry biomass (g)"
            , "HI"  = "Harvest index (HI)" 
            , "SLA" = "Specific leaf area (cm^2^/g)"
            , "RCC" = "Relative chlorophyll content (RCC)"
            , "WUE_B" = "Biomass water use efficiency (g/L)"
            , "WUE_T" = "Tuber water use efficiency (g/L)"
            )

4 Import data

fb <- gs %>% 
  range_read("fb") %>% 
  rename_with("tolower") %>%
  rename_with(~gsub("\\s+|\\.", "_", .)) %>% 
  mutate(treat = case_when(treat %in% "wellwater" ~ "WW"
                           , treat %in% "drydown" ~ "WD")) %>% 
  select(block, treat, genotype,
         spad_29 = spad_29dap, 
         spad_59 = spad_59dap, 
         spad_76 = spad_76dap, 
         spad_83 = spad_83dap,
         hgt = hgt_86dap,
         rwc = rwc_84dap,
         lop = op_84dap,
         ldw = leafdw,
         sdw = stemdw,
         rdw = rootdw,
         tdw = tubdw,
         ntub,
         trs = ttrns,
         lfa = la,
         rtl = rlg
         ) %>% 
  mutate(tdb = (ldw+sdw+rdw+tdw),
         hi = tdw/(ldw+sdw+rdw+tdw),
         sla = lfa/ldw,
         rcc = (spad_83/lfa),
         wue_b = (ldw+sdw+rdw+tdw)/trs,
         wue_t = tdw/trs
         ) %>% 
  mutate(gnt = recode(genotype, !!!geno), .after = genotype) %>% 
  select(genotype, gnt, treat, block, everything()) %>% 
  mutate(across(where(is.character), as.factor)) %>% 
  rename_with(., toupper, where(is.numeric)) 

# str(fb)

fb %>% web_table()

5 Abbreviations

gs %>% 
  range_read("var") %>% 
  filter(include == "x") %>%
  select(Variable, Abbreviation) %>%
  kable()

Variable	Abbreviation
Soil Plant Analysis Development	SPAD
Plant height	HGT
Relative water content	RWC
Leaf osmotic potential	LOP
Leaf dry weight	LDW
Stem dry weight	SDW
Root dry weight	RDW
Tuber dry weight	TDW
Tuber number	NTUB
Root length	RTL
Total transpiration	TRS
Leaf area	LFA
Total dry biomass	TDB
Harvest index	HI
Specific leaf area	SLA
Water use efficiency	WUE
Tuber water use efficiency	TWUE

6 Table 1

tab <- gs %>% 
  range_read("gnt") %>% 
  dplyr::select(Number, Genotypes, Adaptability, "Growing period" = GPL, "Heat tolerance" = Heat, "Dry matter (%)")
  
# tab %>% sheet_write(gs, data = .,  sheet = "tab1")

tab %>% web_table()

7 Table 2

tab <- fb %>%
  select(!c(block, gnt, genotype)) %>% 
  tbl_summary(
    by = treat, 
    statistic = list(all_continuous() ~ "{mean} ± {sd}"),
    missing = "no") %>% 
  add_p() %>% 
  bold_labels() %>% 
  as_tibble() %>% 
  mutate(`**Characteristic**` = recode(
    `**Characteristic**`
    , "__SPAD_29__" = "Chlorophyll concentration (SPAD) at 29 dap"
    , "__SPAD_59__" = "Chlorophyll concentration (SPAD) at 59 dap"
    , "__SPAD_76__" = "Chlorophyll concentration (SPAD) at 76 dap"
    , "__SPAD_83__" = "Chlorophyll concentration (SPAD) at 83 dap"
    , "__HGT__" = "Plant height (cm)"
    , "__RWC__" = "Relative water content (%)"
    , "__LOP__" = "Leaf osmotic potential (MPa)"
    , "__LDW__" = "Leaf dry weight (g)"
    , "__SDW__" = "Stem dry weight (g)"
    , "__RDW__" = "Root dry weight (g)"
    , "__TDW__" = "Tuber dry weight (g)"
    , "__NTUB__" = "Tuber number (N°)" 
    , "__TRS__" = "Total transpiration (mL)"
    , "__LFA__" = "Leaf area (cm^2^)"
    , "__RTL__" = "Root length (cm)"
    , "__TDB__" = "Total dry biomass (g)"
    , "__HI__"  = "Harvest index (HI)" 
    , "__SLA__" = "Specific leaf area (cm^2^/g)"
    , "__RCC__" = "Relative chlorophyll content (RCC)"
    , "__WUE_B__" = "Biomass water use efficiency (g/L)"
    , "__WUE_T__" = "Tuber water use efficiency (g/L)"
    )) %>% 
  rename(Variable = `**Characteristic**`
         , `Water deficit` = `**WD**, N = 75`
         , `Well-Watered`= `**WW**, N = 75`
         , `p-value` = `**p-value**`) 

# tab %>% sheet_write(gs, data = .,  sheet = "tab2")

tab %>% web_table()

8 Figure 1

fts <- gs %>% 
  range_read("ftsw") %>%
  filter(Treatment != "preharvest") %>%
  tidyr::gather(key = day, value = fts, -ID, -Genotype, -Treatment)

model <- fts %>% 
  inti::mean_comparison(response = "fts"
                        , model_factors = "Treatment*day"
                        , c("Treatment", "day"))

plt1 <- model$comparison %>% 
  mutate(Treatment = case_when(
    Treatment == "drydown" ~ "Water Deficit (WD)"
    , Treatment == "wellwater" ~ "Well Watered (WW)"
  )) %>% 
  mutate(across(day, as.numeric)) %>% 
  plot_smr(
  type = "line", color = F
  , x = "day"
  , y = "fts"
  , group = "Treatment"
  , ylab = "Fraction of transpirable soil water"
  , xlab =  "Days"
  , glab = "Treatment"
  , legend = "none"
  , ylimits = c(0, 1.08, 0.1)
  , error =  "ste"
  ) +
  scale_x_continuous(breaks = c(0:100)*2
                     , limits = c(34, 88))

trns <- gs %>% 
  range_read("trans") %>%
  filter(Treatment != "preharvest") %>%
  tidyr::gather(key = day, value = trans, -ID, -Genotype, -Treatment) %>%
  filter(day != "TOTAL") %>%
  drop_na()


model <- trns %>% 
  inti::mean_comparison(response = "trans"
                        , model_factors = "Treatment*day"
                        , c("Treatment", "day"))

plt2 <- model$comparison %>% 
    mutate(Treatment = case_when(
    Treatment == "drydown" ~ "Water Deficit (WD)"
    , Treatment == "wellwater" ~ "Well Watered (WW)"
  )) %>% 
  mutate(across(day, as.numeric)) %>% 
  plot_smr(type = "line", color = F
            , x = "day"
            , y = "trans"
            , group = "Treatment"
            , ylab = "Transpiration (mL)"
            , xlab =  "Days"
            , glab = "Irrigation"
            , ylimits = c(0, 700, 100)
            , error = "ste"
            ) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

plot <- ggdraw(xlim = c(0, 0.5), ylim = c(0,0.9)) +
  draw_plot(plt2, width = 0.5, height = 0.43 , x = 0.0, y = 0.47 ) +
  draw_plot(plt1, width = 0.495, height = 0.5, x = 0.005, y = 0.0 ) +
          draw_plot_label(
            label = c("a", "b"),
            x = c(0.465, 0.465),
            y = c(0.79, 0.49))

plot %>% 
  ggsave2(plot = .
          , "files/Fig1.pdf", width = 18, height = 15, units = "cm")

"files/Fig1.pdf" %>% include_graphics()

9 Figure 2

# tdw

model <- fb %>% 
  inti::mean_comparison(response = "TDW"
                        , model_factors = "block + gnt*treat"
                        , c("gnt","treat"))

tdw <- model$comparison %>% 
  plot_smr(type = "bar"
          , color = F
          , x = "gnt"
          , xlab = ""
          , y = "TDW"
          , ylab = "Tuber dry weight (g)"
          , group = "treat"
          , glab = "Treatment"
          , legend = "none"
          , error = "ste"
          , ylimits = c(0, 100, 20)
          )

# rcc

model <- fb %>% 
  inti::mean_comparison(response = "RCC"
                        , model_factors = "block + gnt*treat"
                        , c("gnt","treat"))

rcc <- model$comparison %>% 
  plot_smr(type = "bar"
        , color = F
        , x = "gnt"
        , xlab =  ""
        , y = "RCC"
        , ylab = "Relative chlorophyll content"
        , ylimits = c(0, 0.1, 0.02)
        , group = "treat"
        , glab = "Treatment"
        , legend = "none"
        , error = "ste"
        ) 

# hi

model <- fb %>% 
  inti::mean_comparison(response = "HI"
                        , model_factors = "block + gnt*treat"
                        , c("gnt","treat"))

hi <- model$comparison %>% 
  plot_smr(type = "bar"
        , color = F
        , x = "gnt"
        , xlab =  ""
        , y = "HI"
        , ylab = "Harvest Index"
        , group = "treat"
        , glab = "Treatment"
        , legend = "none"
        , ylimits = c(0, 1, 0.2)
        , error = "ste"
        ) 

# wue_t

model <- fb %>% 
  inti::mean_comparison(response = "WUE_T"
                        , model_factors = "block + gnt*treat"
                        , c("gnt","treat"))

wue_t <- model$comparison %>% 
  plot_smr(type = "bar", color = F
        , x = "gnt"
        , xlab =  ""
        , y = "WUE_T"
        , ylab = "Tuber water use efficiency (g/L)"
        , group = "treat"
        , glab = "Treatment"
        , legend = "none"
        , ylimits = c(0, 10, 2)
        , error = "ste"
        ) 

# spad_29

model <- fb %>% 
  inti::mean_comparison(response = "SPAD_29"
                        , model_factors = "block + gnt*treat"
                        , c("gnt","treat"))

spad_29 <- model$comparison %>% 
    mutate(treat  = case_when(
    treat  == "WD" ~ "Water Deficit (WD)"
    , treat  == "WW" ~ "Well Watered (WW)"
    )) %>%
  plot_smr(type = "line"
        , x = "gnt"
        , xlab =  ""
        , y = "SPAD_29"
        , ylab = "SPAD at 29 dap"
        , group = "treat"
        , glab = "Treatment"
        , legend = "top"
        , ylimits = c(30, 70, 5)
        , error = "ste"
        , color = F
        )

# spad_83

model <- fb %>% 
  inti::mean_comparison(response = "SPAD_83"
                        , model_factors = "block + gnt*treat"
                        , c("gnt","treat"))

spad_83 <- model$comparison %>% 
  plot_smr(type = "line", color = F
        , x = "gnt"
        , xlab =  ""
        , y = "SPAD_83"
        , ylab = "SPAD at 83 dap"
        , group = "treat"
        , glab = "Treatment"
        , legend= "none"
        , ylimits = c(30, 70, 5)
        , error = "ste"
        ) 

plot <- ggdraw(xlim = c(0.0, 1), ylim = c(0.0, 1)) +
  draw_plot(spad_29 , width = 0.49, height = 0.35, x = 0.0, y = 0.6) +
  draw_plot(spad_83, width = 0.49, height = 0.3, x = 0.5, y = 0.6) +
  draw_plot(rcc,  width = 0.49, height = 0.3, x = 0.0, y = 0.3) +
  draw_plot(tdw, width = 0.49, height = 0.3, x = 0.5, y = 0.3) +
  draw_plot(hi , width = 0.49, height = 0.3, x = 0.0, y = 0.0) +
  draw_plot(wue_t, width = 0.49, height = 0.3, x = 0.5, y = 0.0) +
          draw_plot_label(
            label = c("a", "b", "c", "d", "e", "f"),
            x = c(0.06, 0.56, 0.06, 0.56, 0.06, 0.56),
            y = c(0.89, 0.89, 0.59, 0.59, 0.29, 0.29))

plot %>% 
  ggsave(plot = .
         , "files/Fig2.pdf"
         , units = "cm", width = 28, height = 25)

"files/Fig2.pdf" %>% include_graphics()

10 Rename variables in PCA

vars <- names(fb)

model <- 5:length(vars) %>% map(function(x) {
  
    fb %>% 
      H2cal(trait = vars[x]
            , gen.name = "genotype"
            , rep.n = 4
            , fixed.model = "0 + (1|block) + genotype"
            , random.model = "1 + (1|block) + (1|genotype)"
            , emmeans = T
            , plot_diag = T
            , outliers.rm = T
            )
})

tabsmr <- 1:length(model) %>% map(function(x) {
      model[[x]][["tabsmr"]] 
      }) %>% 
  bind_rows()

varnames <- tabsmr %>% 
  select(trait, repeatability) %>% 
  mutate(name = paste0("(",round(repeatability, 2), ") ", trait)) %>% 
  select(!repeatability) %>% 
  deframe()

11 Figure 3

fig <- magick::image_read("files/potatos.png") %>% 
  grid::rasterGrob() %>%
  ggsave("files/Fig3.pdf", plot =  ., width = 13, height = 13)

"files/Fig3.pdf" %>% include_graphics()

12 Figure 4

mvdata <- fb %>% 
  group_by(genotype, treat) %>% 
  summarise(across(where(is.numeric), ~mean(.x, na.rm = T))) %>% 
  ungroup() %>% 
  unite(trait, c(genotype, treat), sep = "-") 

plot <- tempfile(fileext = ".pdf")
pdf(file = plot, width = 20, height = 20)
mvdata %>% 
  select(where(is.numeric)) %>% 
  select(!ends_with(c("29", "59", "76"))) %>% 
  pairs.panels(x = .
               , hist.col="red"
               , pch = 21
               , stars = TRUE
               , scale = FALSE
               , lm = TRUE
               ) 
graphics.off()

fig <- magick::image_read_pdf(plot) %>% 
  grid::rasterGrob() %>%
  ggsave2("files/Fig4.pdf", plot =  ., width = 20, height = 20, units = "cm")

"files/Fig4.pdf" %>% include_graphics()

13 Figure 5

mv <- mvdata %>% 
  pivot_longer(!trait) %>% 
  mutate(name = recode(name, !!!varnames)) %>% 
  pivot_wider() %>% 
  column_to_rownames("trait") %>% 
  PCA(scale.unit = T, graph = F)

hcpc <- HCPC(mv, graph = F)

outfile <- tempfile(fileext = ".pdf")
pdf(outfile, width = 9, height = 9)
plot.HCPC(x = hcpc
          , choice = "map"
              , legend = list(x = "topright"
                              , cex = 0.6
                              , inset = 0.001
                              , box.lty=0
                              )
              , draw.tree = F
          )
graphics.off()

tree <- magick::image_read_pdf(outfile) %>% 
  grid::rasterGrob()

pca <- mv %>% 
  plot.PCA(choix = "var", autoLab = "auto", cex = 0.6)

plot <- plot_grid(pca, tree
                  , nrow = 1, labels = "auto")
plot %>% 
  ggsave2(plot = .
          , "files/Fig5.pdf"
          , units = "cm", width = 30, height = 15)

"files/Fig5.pdf" %>% include_graphics()

14 Supplementary information

15 Supplementary Table 1

h2plot <- tabsmr %>% 
  mutate(trait = recode(trait, !!!traits)) %>% 
  select(trait, mean, std, min, max, V.g, V.e, repeatability) %>% 
  mutate(across(where(is.numeric), ~ round(.x, digits = 2))) %>% 
  mutate(trait = recode(trait, !!!traits)) %>% 
  rename(Trait = trait, H2 = repeatability)
  
h2plot  %>% web_table()

# h2plot %>% sheet_write(gs, data = .,  sheet = "tabS1")

16 Supplementary Figure 1

var <- get_pca_var(mv)

tmp <- tempfile(fileext = ".png")
ppi <- 300
png(tmp, width=3.8*ppi, height=10*ppi, res=ppi)
corrplot(var$cor, 
         method="number",
         tl.col="black", 
         tl.srt=45,)
graphics.off()

pt1 <- png::readPNG(tmp) %>%
  grid::rasterGrob(interpolate = TRUE)

pt2 <- fviz_eig(mv, 
                addlabels=TRUE,
                hjust = 0.05,
                barfill="white",
                barcolor ="darkblue",
                linecolor ="red") + 
  ylim(0, 50) + 
  labs(
    title = "PCA - percentage of explained variances",
    y = "Variance (%)") +
  theme_minimal()

pt3 <- fviz_contrib(mv,
                     choice = "var", 
                     axes = 1, 
                     top = 10,
                     fill="white",
                     color ="darkblue",
                     sort.val = "desc") +
  ylim(0, 18) + 
  labs(title = "Dim 1 - variables contribution") 

pt4 <- fviz_contrib(mv,
                     choice = "var", 
                     axes = 2, 
                     top = 10,
                     fill="white",
                     color ="darkblue",
                     sort.val = "desc") +
  ylim(0, 18) + 
  labs(title = "Dim 2 - variables contribution") 

plot <- ggdraw(xlim = c(0.0, 1.0), ylim = c(0, 1.0))+
  draw_plot(pt1,  width = 0.4, height = 0.99, x = 0.6, y = 0.0) +  
  draw_plot(pt2,  width = 0.6, height = 0.34, x = 0.03, y = 0.66) +
  draw_plot(pt3, width = 0.6, height = 0.34, x = 0.03, y = 0.33) + 
  draw_plot(pt4, width = 0.6, height = 0.34, x = 0.03, y = 0.0) +
        draw_plot_label(
    label = c("a", "b", "c", "d"),
    x = c(0.005, 0.005, 0.005, 0.65),
    y = c(0.999, 0.67, 0.34, 0.999))

ggsave2(plot = plot, "files/FigS1.pdf", height = 25, width = 30, units = "cm")

"files/FigS1.pdf" %>% include_graphics()

17 Convert figures to pdf to png

figs <- "files/"  %>%  list.files(pattern = "pdf", full.names = T)

convert <- 1:length(figs) %>% map(\(x) {
  figname <- figs[x] %>% basename() %>% gsub("pdf", "png", .)
  
  figs[x] %>% 
    image_read_pdf() %>% 
    image_write(file.path("files", figname))
})